import holoviews as hv
hv.extension('bokeh')
hv.opts.defaults(hv.opts.Curve(width=500),
hv.opts.Image(width=500, colorbar=True, cmap='Viridis'))
import numpy as np
import scipy.signal
import scipy.fft
from IPython.display import Audio
9. Diseño de sistemas y filtros IIR¶
Un filtro FIR de buena calidad puede requerir una gran cantidad de coeficientes
Es posible implementar filtros más eficientes usando recursividad. Esta es la base de los filtros de respuesta al impulso infinita o IIR que veremos en esta lección
9.1. Definición de un sistema IIR¶
Generalizando el sistema FIR para incluir versiones pasadas de la salida y asumiendo \(a[0] = 1\) llegamos a
es decir dos convoluciones discretas que definen una ecuación de diferencias
Este tipo de sistema se conoce como
sistema infinite impulse response (IIR)
sistema auto-regresive moving average (ARMA)
autoregresivo de orden M: incluye valores pasados de la salida
media movil de orden L+1: pondera el valor presente y pasados de la entrada
Podemos ver el sistema IIR como una generalización del sistema FIR. El caso particular del sistema FIR se recupera si
\(a[m] = 0\) para \(m=[1, \ldots, M]\)
9.1.1. Respuesta en frecuencia del sistema IIR¶
Aplicando la transformada de Fourier convertimos las convoluciones en multiplicaciones y encontramos la respuesta en frecuencia como
que existe siempre que \(A[k] \neq 0\).
La respuesta en frecuencia también suele expresarse como
donde
\(K\) se llama ganancia
las raices del polinomio del numerador \(\alpha\) se llaman conjuntamente ceros
las raices del polinomio del denominador \(\beta\) se llaman conjuntamente polos
9.1.2. Ejemplo de respuesta al impulso de un sistema IIR¶
Consideremos el siguiente sistema IIR
Los coeficientes del sistema son
\(a[0] = 1\), \(a[1] = -\gamma\) y \(b[0] = (1-\gamma)\)
Es decir que es AR de orden 1 y MA de orden 1
¿Cúal es su respuesta al impulso? Asumiendo \(y[n]=0, n<0\), tenemos que
¿Cómo cambia la respuesta al impulso con distintos valores de \(\gamma\)? ¿Qué pasa si \(\gamma \geq 1\)?
Respondamos estas preguntas visualizando la respuesta al impulso de este sistema con la función scipy.signal.dimpulse
# Valores de gamma que probaremos:
gamma = [-1.5, -1, -0.5, 0.5, 1., 1.5]
p = []
for g in gamma:
t, y = scipy.signal.dimpulse(([1-g, 0], [1,-g], 1), x0=0, n=30)
p.append(hv.Curve((t, y[0][:, 0]), label=f"gamma={g}"))
hv.Layout(p).cols(3).opts(hv.opts.Curve(width=250, height=200, axiswise=True))
/home/phuijse/.conda/envs/info320/lib/python3.9/site-packages/scipy/signal/filter_design.py:1630: BadCoefficients: Badly conditioned filter coefficients (numerator): the results may be meaningless
warnings.warn("Badly conditioned filter coefficients (numerator): the "
De las figuras podemos ver que:
Para \(\gamma < 0\) (primera fila) los coeficientes del sistema son alternantes en signo
Para \(|\gamma| < 1\) los coeficientes del sistema tienden a cero
Para \(|\gamma| > 1\) los coeficientes del sistema divergen y tienen a infinito
Advertencia
A diferencia de un sistema FIR, el sistema IIR puede tener configuraciones inestables en que los coeficientes crecen o decrecen infinitamente
Por otro lado consideremos el sistema anterior y asumamos que \(|\gamma|<1\), desenrollando tenemos que
Nota
Con un sistema IIR de pocos coeficientes podemos representar un sistema FIR considerablemente más grande
En el ejemplo anterior, si escogemos \(\gamma\) tal que \(\gamma^{20 }\approx 0\) entonces aproximamos un sistema FIR de orden 20 con tan sólo 3 coeficientes
9.1.3. Ejemplo de respuesta en frecuencia de un sistema IIR¶
Para el sistema del ejemplo anterior su respuesta en frecuencia es
que en notación de polos y ceros se escribe como
es decir que tiene un cero en \(0\), un polo en \(\gamma\) y una ganancia de \((1-\gamma)\)
Para entender mejor este sistema estudiemos la magnitud de \(|H[k]|\) para \(\gamma < 1\)
¿Cómo se ve \(|H[k]|\)? ¿Qué función cumple este sistema?
k = np.arange(-24, 25)/50
Hk = lambda gamma, k : (1-gamma)/np.sqrt(1 - 2*gamma*np.cos(2.0*np.pi*k) + gamma**2)
p = []
for gamma in [0.25, 0.5, 0.75]:
p.append(hv.Curve((k, Hk(gamma, k)), 'Frecuencia', 'Respuesta', label=f'gamma={gamma}'))
hv.Overlay(p)
Nota
Este sistema atenua las frecuencias altas, es decir que actua como un filtro pasa bajos
9.2. Diseño de filtros IIR simples¶
Los filtros IIR más simples son los de un polo y un cero, es decir filtros de primer orden
donde podemos reconocer
\(b[0]=K\)
\(\beta = - b[1] \cdot K\)
\(\alpha=-a[1]\)
Definimos la frecuencia de corte \(f_c\) como aquella frecuencia en la que el filtro alcanza una atenuación de 0.7 (-3 dB). Haciendo la equivalencia con el ejemplo anterior tenemos que \(\gamma = e^{-2\pi f_c}\)
9.2.1. Receta para un filtro pasa bajo IIR con frecuencia de corte \(f_c\)¶
Asignamos
\(b[0] = 1 - e^{-2\pi f_c}\)
\(b[1] = 0\)
\(a[1] = -e^{-2\pi f_c}\)
Lo que resulta en la siguiente respuesta en frecuencia
Es decir un cero en \(0\), un polo en \(e^{-2\pi f_c}\) y ganancia \(1-e^{-2\pi f_c}\)
9.2.2. Receta para un filtro pasa alto IIR con frecuencia de corte \(f_c\)¶
Asignamos
\(b[0] = (1 + e^{-2\pi f_c})/2\)
\(b[1] = -(1 + e^{-2\pi f_c})/2\)
\(a[1] = -e^{-2\pi f_c}\)
Lo que resulta en la siguiente respuesta en frecuencia
Es decir un cero en \(1\), un polo en \(e^{-2\pi f_c}\) y ganancia \(\frac{1+e^{-2\pi f_c}}{2}\)
9.2.3. Aplicar un filtro a una señal con scipy¶
Para filtrar una señal unidimensional con un filtro IIR (sin variar la fase de la señal) podemos utilizar la función
scipy.signal.filtfilt(b, # Coeficientes del numerador
a, # Coeficientes del denominador
x, # Señal a filtrar
...
)
Los siguientes ejemplos muestran un señal de tipo pulso rectangular filtrada con sistemas IIR de primer orden pasa bajo y pasa-alto diseñados con las recetas mostradas anteriormente
n = np.arange(0, 500)
x = 0.5 + 0.5*scipy.signal.square((n)/(2.*np.pi*5), duty=0.3)
def iir_low_pass(signal, fc):
gamma = np.exp(-2*np.pi*(fc))
b, a = [(1-gamma), 0], [1, -gamma]
return scipy.signal.filtfilt(b, a, signal)
y = {}
for fc in [0.05, 0.02, 0.01]:
y[fc] = iir_low_pass(x, fc)
px = hv.Curve((n, x))
py = []
for fc, y_ in y.items():
py.append(hv.Curve((n, y_), label=f'fc={fc}'))
hv.Layout([px, hv.Overlay(py)]).cols(1).opts(hv.opts.Curve(height=200))
def iir_high_pass(signal, fc):
gamma = np.exp(-2*np.pi*(fc))
b, a = [(1+gamma)/2, -(1+gamma)/2], [1, -gamma]
return scipy.signal.filtfilt(b, a, signal)
y = {}
for fc in [0.01, 0.02, 0.05]:
y[fc] = iir_high_pass(x, fc)
px = hv.Curve((n, x))
py = []
for fc, y_ in y.items():
py.append(hv.Curve((n, y_), label=f'fc={fc}'))
hv.Layout([px, hv.Overlay(py)]).cols(1).opts(hv.opts.Curve(height=200))
Nota
El filtro pasa-bajos suaviza los cambios de los pulsos rectangulares. El filtro pasa-altos elimina las zonas constantes y resalta los cambios de la señal.
9.3. Diseño de filtros IIR de segundo orden¶
Los filtros IIR de segundo orden o biquad tienen dos polos y dos ceros.
Su respuesta en frecuencia es
donde \(W_N = e^{-j \frac{2 \pi}{N}}\) y la relación entreo coeficientes y polos/ceros es:
Con arquitecturas de segundo orden se pueden crear filtros pasabanda y rechaza banda
9.4. Diseño de filtros IIR de orden mayor¶
Para crear los coeficientes de filtro IIR de orden mayor podemos usar la función
scipy.signal.iirfilter(N, # Orden del filtro
Wn, # Frecuencias de corte (normalizadas en [0,1])
fs, # Frecuencia de muestreo
btype='bandpass', # Tipo de filtro: 'bandpass', 'lowpass', 'highpass', 'bandstop'
ftype='butter', # Familia del filtro: 'butter', 'ellip', 'cheby1', 'cheby2', 'bessel'
output='ba', # Retornar coeficientes
...
)
El filtro Butterworth es óptimo en el sentido de tener la banda de paso lo más plana posible.
Otros filtros se diseñaron con otras consideraciones.
Los filtros IIR digitales están basados en los filtros IIR analógicos.
Observe como al aumentar el orden el filtro pasabajo IIR comienza a cortar de forma más abrupta
Hk = {}
for order in [1, 2, 5, 20]:
b, a = scipy.signal.iirfilter(N=order, Wn=0.2, fs=1,
ftype='butter', btype='lowpass', output='ba')
freq, response = scipy.signal.freqz(b, a, fs=1)
Hk[order] = np.abs(response)
p = []
for order, response in Hk.items():
p.append(hv.Curve((freq, response), 'Frecuencia', 'Respuesta', label=f'orden={order}'))
hv.Overlay(p)
9.5. Comparación de la respuesta en frecuencia de filtros FIR e IIR del orden equivalente¶
Comparemos la respuesta en frecuencia de un filtro IIR y otro FIR ambos pasa-bajo con 20 coeficientes
Fs = 1
fc = 0.25
h = scipy.signal.firwin(numtaps=20, cutoff=fc, pass_zero=True, window='hann', fs=Fs)
b, a = scipy.signal.iirfilter(N=9, Wn=fc, fs=Fs, ftype='butter', btype='lowpass')
display(len(h), len(b)+len(a))
freq_fir, response_fir = scipy.signal.freqz(h, 1, fs=Fs)
freq_iir, response_iir = scipy.signal.freqz(b, a, fs=Fs)
20
20
p1 = hv.Curve((freq_fir, np.abs(response_fir)), 'Frecuencia', 'Respuesta', label='FIR')
p2 = hv.Curve((freq_iir, np.abs(response_iir)), 'Frecuencia', 'Respuesta', label='IIR')
hv.Overlay([p1, p2])*hv.VLine(fc).opts(color='k', alpha=0.5)
La linea negra marca la ubicación de la frecuencia de corte
Nota
El filtro IIR es mucho más abrupto, es decir filtra mejor, que el filtro FIR equivalente
Una desventaja del filtro IIR es que por definición introduce una desfase no constante en la señal de salida
freq_fir, delay_fir = scipy.signal.group_delay(system=(h, 1), fs=Fs)
freq_iir, delay_iir = scipy.signal.group_delay(system=(b, a), fs=Fs)
/home/phuijse/.conda/envs/info320/lib/python3.9/site-packages/scipy/signal/filter_design.py:688: UserWarning: The group delay is singular at frequencies [3.105, 3.111, 3.117, 3.123, 3.129, 3.135], setting to 0
warnings.warn(
p1 = hv.Curve((freq_fir, delay_fir), 'Frecuencia', 'Desfase', label='FIR')
p2 = hv.Curve((freq_iir, delay_iir), 'Frecuencia', 'Desfase', label='IIR')
hv.Overlay([p1, p2])*hv.VLine(fc).opts(color='k', alpha=0.5)
¿Cómo se ve una señal filtrada donde se preserva la fase versus una donde no se preserva la fase?
Consideremos la señal rectangular anterior y apliquemos un filtro pasa-bajo IIR de orden 1
Esta vez compararemos el filtro con la función scipy.signal.lfilter y la función scipy.signal.filtfilt. La primera no preserva la fase mientras que la segunda si lo hace
Fs = 1
fc = 0.01
n = np.arange(0, 500)
x = 0.5 + 0.5*scipy.signal.square((n)/(2.*np.pi*5), duty=0.3)
b, a = scipy.signal.iirfilter(N=1, Wn=fc, fs=Fs, ftype='butter', btype='lowpass')
# No se preserva la fase
y_lfilter = scipy.signal.lfilter(b, a, x)
# Se preserva la fase
y_filtfilt = scipy.signal.filtfilt(b, a, x)
px = hv.Curve((n, x), 'Tiempo', 'Entrada')
py = []
py.append(hv.Curve((n, y_filtfilt), 'Tiempo', 'Salida', label=f'Fase constante'))
py.append(hv.Curve((n, y_lfilter), 'Tiempo', 'Salida', label=f'Fase no constante'))
hv.Layout([px, hv.Overlay(py)]).cols(1).opts(hv.opts.Curve(height=200))
Nota
En el caso donde no se preserva la fase podemos notar que la señal de salida está desplazada con respecto a la original. Además los cambios tienen una transición asimétrica
La función scipy.signal.filtfilt “arregla” el problema del desfase filtrando la señal dos veces. La primera vez se filtra hacia adelante en el tiempo y la segunda vez hacia atrás. Por ende no se puede aplicar en un escenario de tipo streaming donde los datos van llegando de forma causal.
En una aplicación causal donde se necesite preservar la fase debemos usar un filtro FIR.
9.6. Apéndice: Efectos de audio con filtros IIR¶
El siguiente ejemplo muestra como implementar el conocido filtro Wah-wah usando un sistema IIR
Este es un filtro pasabanda modulado con ancho de pasada fijo \(f_b\) [Hz] y una frecuencia central variable \(f_c\) [Hz], donde La frecuencia central se modula con una onda lenta
Se modela como el siguiente sistema IIR
donde
y
Veamos como modifica este filtro una señal de audio
import librosa
data, fs = librosa.load("../../data/DPSAU.ogg")
Audio(data, rate=fs)
data_wah = []
zi = np.zeros(shape=(2,))
# Parámetros fijos del filtro
fb, Nw = 200, 5
c = (np.tan(np.pi*fb/fs) - 1.)/(np.tan(2*np.pi*fb/fs) +1)
# Filtramos una ventana de la señal moviendo lentamente fc
for k in range(len(data)//Nw):
# Cálculo de la frecuencia central
fc = 500 + 2000*(np.cos(2.0*np.pi*k*30./fs) +1)/2
d = -np.cos(2*np.pi*fc/fs)
# Coeficientes del filtro
b, a = [(1+c), 0, -(1+c)], [1, d*(1-c), -c]
# Filtramos, usando el filtrado anterior como borde (zi)
data2, zi = scipy.signal.lfilter(b, a, data[k*Nw:(k+1)*Nw], zi=zi)
# Guardamos
data_wah.append(data2)
Audio(np.hstack(data_wah), rate=int(fs))
Si quieres profundizar en el tema de los filtros IIR aplicados a efectos de audio recomiendo: https://www.ee.columbia.edu/~ronw/adst-spring2010/lectures/lecture2.pdf